function [Param, Covar, Resid, Info] = ecmlsrmle(Data, Design, ...
   MaxIter, TolParam, TolObj, Param0, Covar0, CovarFormat)
%ECMLSRMLE Least-squares regression (with missing data).
% Estimate a NUMPARAMS x 1 column vector of model parameters Param and a NUMSERIES x NUMSERIES
% residual covariance matrix Covar for a least-squares regression model with missing data, where the
% model has the equivalent form
%		Data(k)  ~  N(Design(k) * Param, Covar)
% for samples k = 1, ... , NUMSAMPLES. If no output arguments, display a plot of convergence (or
% non-convergence) of algorithm to a maximum likelihood estimate.
%
%		[Param, Covar] = ecmlsrmle(Data, Design);
%		[Param, Covar, Resid, Info] = ecmlsrmle(Data, Design, ...
%			MaxIter, TolParam, TolObj, Param0, Covar0, CovarFormat);
%
%		ecmlsrmle(Data, Design);
%		ecmlsrmle(Data, Design, MaxIter, TolParam, TolObj, Param0, Covar0, CovarFormat);
%
% Inputs:
%	Data - NUMSAMPLES x NUMSERIES matrix with NUMSAMPLES samples of a NUMSERIES-dimensional random
%		vector. Missing values are represented as NaNs. Only samples that are entirely NaNs are
%		ignored (to ignore samples with at least one NaN, use MVNRMLE).
%	Design - Either a matrix or a cell-array to handle two distinct model structures. First, if
%		NUMSERIES = 1, Design can be a NUMSAMPLES x NUMPARAMS matrix with known values. This is the
%		"standard" form for regression on a single data series. Alternatively, for any
%		NUMSERIES >= 1, Design can be a cell array of either 1 or NUMSAMPLES cells, where each cell
%		contains a NUMSERIES x NUMPARAMS matrix of known values. If Design has a single cell, then
%		it is assumed to be the same Design matrix for each sample. Otherwise, Design must contain
%		individual Design matrices for each sample.
%
% Optional Inputs:
%	MaxIter - Maximum number of iterations for the estimation algorithm. Default value is 100.
%	TolParam - Convergence tolerance for estimation algorithm based on changes in model parameter
%		estimates. Default value is sqrt(eps) which is about 1.0e-8 for double precision. The
%		convergence test for changes in model parameters is
%			||Param(k) - Param(k-1)|| < TolParam * (1 + ||Param(k)||)
%		for iteration k = 2, 3, ... . Convergence is assumed when both the TolParam and TolObj
%		conditions are satisfied. If both TolParam <= 0 and TolObj <= 0, then do maximum number of
%		iterations (MaxIter) regardless of the results of convergence tests.
%	TolObj - Convergence tolerance for estimation algorithm based on changes in the objective
%		function. Default value is eps^(3/4) which is about 1.0e-12 for double precision. The
%		convergence test for changes in the objective function is
%			|Obj(k) - Obj(k-1)| < TolObj * (1 + |Obj(k)|)
%		for iteration k = 2, 3, ... . Convergence is assumed when both the TolParam and TolObj
%		conditions are satisfied. If both TolParam <= 0 and TolObj <= 0, then do maximum number of
%		iterations (MaxIter) regardless of the results of convergence tests.
%	Param0 - NUMPARAMS x 1 column vector that contains a user-supplied initial estimate for the
%		parameters of the regression model. Default is a zero vector.
%	Covar0 - NUMSERIES x NUMSERIES matrix that contains a user-supplied initial or known estimate
%		for the covariance matrix of the residuals of the regression. Default is an identity matrix.
%		If doing covariance-weighted least-squares (CWLS), this matrix yields weights for each
%		series in the regression in addition to its function as an initial guess of the residual
%		covariance for the ECM algorithm.
%	CovarFormat - String that specifies the format for the covariance matrix. The choices are:
%		'full' - Default method. Estimate the full covariance matrix.
%		'diagonal' - Estimate the covariance matrix as a diagonal matrix.
%
% Outputs:
%	Param - NUMPARAMS x 1 column vector of estimates for the model parameters of the regression
%		model.
%	Covar - NUMSERIES x NUMSERIES matrix of estimates for the covariance of the residuals of the
%		regression (note that with least-squares, this estimate may not be a maximum likelihood
%		estimate except under special circumstances).
%	Resid - NUMSAMPLES x NUMSERIES matrix of residuals from the regression. For any missing values
%		in Data, the corresponding residual is the difference between the conditionally-imputed
%		value for Data and the model, i.e., the imputed residual. WARNING: The covariance estimate
%		Covar cannot be derived from the residuals.
%	Info - A structure that contains additional information from the regression. The structure has
%		the following fields:
%		Info.Obj - A variable-extent column vector with no more than MAXITER elements that contains
%			each value of the objective function at each iteration of the estimation algorithm. The
%			last value in this vector, which is Obj(end), is the terminal estimate of the objective
%			function. If doing least-squares, the objective function is the sum of the squares of
%			the residuals.
%		Info.PrevParam - NUMPARAMS x 1 column vector of estimates for the model parameters from the
%			iteration just prior to the terminal iteration.
%		Info.PrevCovar - NUMSERIES x NUMSERIES matrix of estimates for the covariance parameters
%			from the iteration just prior to the terminal iteration.
%
% Notes:
%	If doing CWLS, Covar0 should usually be a diagonal matrix. Series with greater influence should
%	have smaller diagonal elements in Covar0 and series with lesser influance shoule have larger
%	diagonal elements. Note that if doing CWLS, Covar0 need not be a diagonal matrix even if
%	CovarFormat = 'diagonal'.
%
%	Although the Design array should not have NaN values, ignored samples due to NaN values in Data
%	are also ignored in the corresponding Design array.
%
%	Note, however, that no NaN values are permitted in Design if it is a 1 x 1 cell array to be used
%	as a single Design matrix for each sample. In addition, a model in this form must have
%	NUMSERIES >= NUMPARAMS and rank(Design{1}) = NUMPARAMS.
%
%	Finally, note that the routine mvnrmle is less strict about the presence of NaN values in the
%	Design array.
%
%	The optional output structure Info with fields Obj, PrevParam, and PrevCovar, is available for
%	diagnostic purposes.
%
% References:
%	[1] Roderick J. A. Little and Donald B. Rubin, Statistical Analysis with Missing Data, 2nd ed.,
%		John Wiley & Sons, Inc., 2002.
%	[2] Xiao-Li Meng and Donald B. Rubin, "Maximum Likelihood Estimation via the ECM Algorithm,"
%		Biometrika, Vol. 80, No. 2, 1993, pp. 267-278.
%	[3] Joe Sexton and Anders Rygh Swensen, "ECM Algorithms that Converge at the Rate of EM,"
%		Biometrika, Vol. 87, No. 3, 2000, pp. 651-662.
%	[4] A. P. Dempster, N.M. Laird, and D. B. Rubin, "Maximum Likelihood from Incomplete Data via
%		the EM Algorithm," Journal of the Royal Statistical Society, Series B, Vol. 39, No. 1, 1977,
%		pp. 1-37.
%
%	See also ECMLSROBJ, ECMMVNRMLE, MVNRMLE

%	Copyright 2005-2007 The MathWorks, Inc.
%	$Revision: 1.1.6.4 $ $Date: 2007/07/03 20:42:08 $

% Step 1 - check arguments

if nargin < 8 || isempty(CovarFormat)
   CovarFormat = 'full';
else
   if ~any(strcmpi(CovarFormat,{'full','diagonal'}))
      error('Finance:ecmlsrmle:InvalidCovarianceFormat', ...
         'Invalid format specified for covariance matrix.');
   end
end
if nargin < 7
   Covar0 = [];
end
if nargin < 6
   Param0 = [];
end
if nargin < 5 || isempty(TolObj)
   TolObj = eps^(3/4);
end
if nargin < 4 || isempty(TolParam)
   TolParam = sqrt(eps);
end
if nargin < 3 || isempty(MaxIter)
   MaxIter = 100;
elseif MaxIter < 1
   MaxIter = 1;
end
if nargin < 2
   error('Finance:ecmlsrmle:MissingInputArg', ...
      'Missing required arguments Data or Design.');
end

[NumSamples, NumSeries, NumParams] = checkmvnrsetup(Data, Design, Param0, Covar0);

if iscell(Design) && (numel(Design) == 1)
   SingleDesign = true;
else
   SingleDesign = false;
end

% Step 2 - observability and ignorability tests

Count = sum(all(isnan(Data),2));
if ((NumSamples - Count) * NumSeries) <= NumParams
   error('Finance:ecmlsrmle:InsufficientData', ...
      'Insufficient data to estimate least-squares regression model.');
end

% Step 3 - initialization

if isempty(Param0)
   Param = zeros(NumParams, 1);
else
   Param = Param0;
end

if isempty(Covar0)						% setup covariance-weighted least-squares
   CovarCWLS = eye(NumSeries, NumSeries);
   Covar = CovarCWLS;
else
   CovarCWLS = Covar0;
   Covar = Covar0;
end

if strcmpi(CovarFormat,'diagonal');		% need to make sure initial covar is diagonal if
   Covar = diag(diag(Covar));			% CovarFormat = 'diagonal' or else premature mle
end

Obj = [];

[CholCovarCWLS, CholState] = chol(CovarCWLS);
if CholState > 0
   warning('Finance:ecmlsrmle:NonPosDefCovar', ...
      'Input covariance matrix is not positive-definite. Will use an identity matrix.');
   CovarCWLS = eye(NumSeries, NumSeries);
   CholCovarCWLS = eye(NumSeries, NumSeries);
end

Resid = nan(NumSamples, NumSeries);

% Step 4 - main loop

for Iter = 1:MaxIter

   if nargout > 3
      Info.PrevParam = Param;
      Info.PrevCovar = Covar;
   end

   LeftSide = zeros(NumParams, NumParams);
   RightSide = zeros(NumParams, 1);
   Z = zeros(NumSeries,1);

   % Step 5 - parameter combined E and CM step

   WarnState = warning('off','MATLAB:nearlySingularMatrix');
   for i = 1:NumSamples
      if iscell(Design)
         if SingleDesign
            Mean = Design{1} * Param;
         else
            Mean = Design{i} * Param;
         end
      else
         Mean = Design(i,:) * Param;
      end

      [mX, mY, CXX, CXY, CYY] = ecmpart(Data(i,:), Mean, Covar);

      if ~isempty(mY)
         if isempty(mX)
            Z = Data(i,:)';
         else
            P = isnan(Data(i,:));
            Q = ~P;
            Y = Data(i,Q)';

            Z(P) = mX + CXY * (CYY \ (Y - mY));
            Z(Q) = Y;
         end

         if iscell(Design)
            if SingleDesign
               A = CholCovarCWLS' \ Design{1};
            else
               A = CholCovarCWLS' \ Design{i};
            end
         else
            A = CholCovarCWLS' \ Design(i,:);
         end
         B = CholCovarCWLS' \ Z;
         LeftSide = LeftSide + (A' * A);
         RightSide = RightSide + (A' * B);
      end
   end
   warning(WarnState);

   Param = LeftSide \ RightSide;

   % Step 6 - covariance E and CM step

   Z = zeros(NumSeries,1);

   Covar0 = Covar;
   Covar = zeros(NumSeries,NumSeries);

   WarnState = warning('off','MATLAB:nearlySingularMatrix');
   Count = 0;
   for i = 1:NumSamples
      if iscell(Design)
         if SingleDesign
            Mean = Design{1} * Param;
         else
            Mean = Design{i} * Param;
         end
      else
         Mean = Design(i,:) * Param;
      end

      [mX, mY, CXX, CXY, CYY] = ecmpart(Data(i,:), Mean, Covar0);
      CovAdj = zeros(NumSeries, NumSeries);

      if ~isempty(mY)
         if isempty(mX)
            Z = Data(i,:)';
         else
            P = isnan(Data(i,:));
            Q = ~P;
            Y = Data(i,Q)';

            Z(P) = mX + CXY * (CYY \ (Y - mY));
            Z(Q) = Y;

            CovAdj(P,P) = CXX - CXY * inv(CYY) * CXY';
         end

         Resid(i,:) = (Z - Mean)';

         Count = Count + 1;
         Covar = Covar + (Resid(i,:)' * Resid(i,:)) + CovAdj;
      end
   end
   warning(WarnState);

   if strcmpi(CovarFormat,'diagonal');
      Covar = diag(diag(Covar));
   end

   Covar = (1.0/Count) .* Covar;

   % Step 7 - evaluate objective and test for convergence

   Objective = ecmlsrobj(Data, Design, Param, CovarCWLS);
   Obj = [Obj; Objective];

   if Iter > 1
      TestObj = Objective - Objective0;
      TestParam = norm(Param - Param0)/sqrt(NumParams);

      EpsObj = TolObj * (1 + abs(Objective));
      EpsParam = TolParam * (1 + norm(Param));

      if ((TestObj >= 0.0) && (TestObj < EpsObj)) && (TestParam < EpsParam)
         break
      end
   end

   Objective0 = Objective;
   Param0 = Param;

   if (Iter == MaxIter) && (MaxIter > 1) && (TolObj > 0.0) && (TolParam > 0.0)
      warning('Finance:ecmlsrmle:EarlyTermination', ...
         'Maximum iterations completed. Convergence criterion not satisfied.');
      break;
   end
end

if nargout > 3
   Info.Obj = Obj;
end

% Step 8 - plot objective function if no output args

if ~nargout
   % Search for existing figure
   findfig = findall(0, 'Type', 'figure', 'tag', 'ecmlsrmle');

   if ~isempty(findfig)
      % Use existing figure
      f = figure(findfig);

   else
      % Create a new figure
      f = figure('Visible', 'Off', 'Tag', 'ecmlsrmle');
   end

   % Plot into current figure
   plot(Obj,'-bo','MarkerFaceColor','b','MarkerSize',3);
   title('\bfProgress of ECM Algorithm in ecmlsrmle');
   xlabel('\bfIteration');
   ylabel('\bfLog-Likelihood');

   % Turn on visibility (if it was off)
   set(f, 'Visible', 'On')
end
